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ABSTRACT 

We investigate the obliquity and spin period of Earth-Moon like systems after 4.5 Gyr of 
tidal evolution with various satellite masses (m s = 0.002577i p - 0.05to p where m p is the 
planet mass) and initial planetary obliquity (e = 0° - 175°) and discuss their relations to 
the habitability of the planet. The satellite initially orbits in the planet's equatorial plane at 
<~ 4 planetary radii and the planet's initial rotation period is 5 h. The other tidal parameters 
are modelled after the Earth and Moon and we keep the satellite on a circular orbit. We find 
three possible outcomes: either i) the system is still evolving, such as our own, ii) the system 
is in the double synchronous state, with the planet's obliquity at either 0° or 180°, or iii) the 
satellite has collided with the planet. The case iii) occurs for initial planetary spins in the 
range Eq ~ 60° - 120°. For other Eq, the satellite survives. The transition between case i) and 
ii) is abrupt and occurs at slightly larger satellite mass (m s ~ 0.02m p ) than the lunar mass. 
For higher masses the system is in the double synchronous state and the final planetary spin 
periods (P p ) are longer than 96 h. We also discuss the habitability of the planet in each case. 
We suggest that cases ii) and iii) are less habitable than case i). Using results from models of 
giant impacts and satellite accretion, we found that the systems that mimic our own i.e., with 
rotation period 12 < P p < 48 h and current planetary obliquity e p < 40° or e p > 140° only 
represent 14% of the possible outcomes. This esimate may only be reliable to within factors 
of a few, depending on how the probability is evaluated. Elser et al. (201 1) conclude that the 
probability of a terrestrial planet having a heavy satellite is 13%. Combining these results 
suggests that the probability of ending up with a system such as our own is of the order of 2%. 

Key words: planets and satellites: general; planets and satellites: dynamical evolution and 
stability; planets and satellites: formation 



1 INTRODUCTION 

With the recent detections of possible terrestrial planets in the 
habitable zones of other stars by the Kepler mission (Lissauer et al., 
2012) a question arises: to what extent could these remote worlds 
support life similar to what we know here on Earth (terrestrial-type 
life)? There are many conditions for a planet to fulfil to support 
terrestrial-type life, which can be controlled by either geological, 
climatological, orbital, or geophysical processes. In this study, 
the first in a series, we focus on the dynamical factors that may 
affect the habitability of a terrestrial exoplanet. We identify several 
key features that we believe are important, if not essential, for 
supporting terrestrial-type life. These are i) a stable climate, ii) 
low to moderate (< 50°) seasonal temperature variations, iii) 
low to moderate diurnal temperature variations, and iv) low to 
moderate spatial variation of temperature over the planet. These 
conditions need to be supplemented with the following: v) regular, 
low-amplitude obliquity oscillations, vi) moderate obliquity to 



induce seasonal variations, vii) moderate rotation rate, and viii) 
small orbital perturbations. We examine each of these below. 

For the current Earth, the first condition requires a constant 
obliquity and a low eccentricity. Earth's life is both land-based 
and water-based and requires a stable climate over millions or 
even billions of years. In addition to continental drift the long-term 
climate on the Earth is largely regulated by the Milankovic 
cycles (Milankovic, 1941). These cycles are related to variations 
in Earth's orbit and their influence on the obliquity. The main 
cycles are related to its equatorial precession (with period 26 kyr, 
corresponding to a frequency ip =-50.5 "/yr, where ip is the angle 
between the vernal equinox and the intersection of the equator 
with the ecliptic), the obliquity variations of the Earth (with period 
41 kyr, corresponding to the frequency ip — S3, where S3 is the 
Earth's nodal eigenfrequency), and the orbital eccentricity (with 
periods 100 kyr and 400 kyr, corresponding to frequencies Q2 + gs 
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and <?2 — gs, where g2 and <?s are the eccentricity eigenfrequencies 
of Venus and Jupiter respectively). The present obliquity variations 
of the Earth are small: the amplitude is just 1.5°, and the low 
amplitude is a result of the presence of the Moon. However, 
even these small oscillations are enough to cause regular ice ages 
through the positive feedback effect (Deser et al., 2000). 

The second condition that is necessary to make the planet 
desirable for terrestrial-type life leans on the following. Spiegel 
et al. (2009) have shown that the polar regions of a planet with 
perpendicular spin are mostly uninhabitable because of the lack 
of seasons: the cold regions will mostly stay cold, possibly even 
below freezing. On Earth the onset of the ice ages tends to favour 
low obliquity for two reasons: the reduction in overall summer 
insolation at high latitudes and the corresponding reduction in 
mean insolation (Huybers & Wunsch, 2005). The cooling would 
result in more ice building up near the polar caps, increasing 
the albedo which then instigates a positive feedback effect. 
However, there are other theories that show that precession and 
eccentricity forcing on the insolation also play a role (e.g. Imbrie 
& Imbrie, 1980; Paillard, 1997; Huybers, 2011) and the debate 
is ongoing. The weak seasonal effects at low obliquity and the 
corresponding constant low temperatures at high latitudes may 
imply that habitability is increased for planets that have a moderate 
obliquity because the associated seasons make the polar regions 
partly habitable by increasing the yearly mean insolation, and 
thus the temperature. However, too high an obliquity may not 
be favourable for habitability either (Williams & Kasting 1997; 
Williams & Pollard, 2003; Spiegel et al., 2009) because the long 
summers at the poles could cause large seasonal temperature 
variations above large continents (Williams & Pollard, 2003) 
and the highest temperatures could exceed 50 C. Additionally 
the long periods of darkness even at mid latitudes could cause 
difficulties for photosynthetic life to survive. Thus it appears that 
the current obliquity favours the development of terrestrial-type 
life rather than more extreme values. 

The third condition, a low to moderate diurnal temperature 
variation, is regulated by a sufficiently slow rotation rate, a thick 
atmosphere and oceans. Diurnal temperature changes are partially 
governed by the relaxation time scale for atmospheric heat losses 
and temperature variations. On the Earth this time scale is 20 
days (Matsuda, 2000) and the diurnal temperature variation is 
approximately 10 K. For Mars the story is very different: the 
relaxation time is comparable to its rotation period and the diurnal 
temperature variation is approximately 60 K (Matsuda, 2000). 
Part of the reason the diurnal variation on Mars is much higher 
than on Earth is because the latter's oceans store heat much more 
effectively than land and Mars' atmosphere is much more tenuous. 
However, a very slow rotation will cause the Hadley cell to cover 
the whole planet (Farrell, 1990) and the diurnal temperature 
variations may be diminished. A fast rotation, on the other hand, 
would decrease the heat transport as argued in the next paragraph. 
By themselves these temperature changes may not pose a problem 
but, coupled with the seasonal variations, they may be difficult for 
land-based life to adapt to. 

The fourth condition, low to moderate spatial variation of 
temperature over the planet, can be the result of the following. Sim- 
ulations of global circulation models on an Earth-like planet with 
different rotation rates have shown that for a fast rotating Earth the 
atmosphere experiences the creation of many small Hadley-like 



cells (Williams, 1988a) and most of the heat transport is caused 
by baroclinic eddies. These eddies decrease the efficiency of heat 
transportation from the warmer regions to the colder regions. 
This possible reduction in heat transportation by the atmosphere 
may decrease the habitability of the planet because the polar (and 
thus colder) regions would probably stay cold (Williams, 1988a). 
However, even though there appears evidence for a more equable 
climate having existed during the Cretaceous and Eocene aeons 
(e.g. Barron, 1989), which may have been caused by latitudinal 
ocean currents rather than the present-day longitudinal ones (Bice 
& Marotzke, 2002), it is not clear to what extent the size of the 
Hadley cells depends on rotation rate or on land mass distribution, 
and how this increased equator to poleward heat flow could have 
been sustained (Barron, 1983). Nevertheless it is probable the area 
of habitability is smaller on a fast-rotating planet than it is on a 
slow-spinning one. Secondly, we have the issue that pertains to 
the effect of tides on ocean mixing. It is well known that tidal 
forcing causes the mixing of stratified layers in the ocean (e.g. 
Egbert & Ray, 2000), though recent work shows that mixing by 
swimming sea creatures could be equally important (Katija & 
Dabiri, 2009). Whatever its source this mixing has profound effects 
on the Earth's climate because it allows for more efficient energy 
exchange between the atmosphere and the ocean, when cooler 
water reaches the surface from the deeper regions. The mixing 
also aids in the transport of water from warmer regions to colder 
ones by currents such as the Gulfstream (Garrett, 2003). Last, the 
mixing brings up important nutrients from the depths of the ocean 
that micro-organisms residing closer to the surface can feed on. 
Thus, it appears that in a simplistic sense ocean mixing could be an 
important ingredient in the sustenance of terrestrial-type life, and 
the habitability of the Earth might be decreased without this effect. 

Regarding the second set of conditions for the existence of 
terrestrial-type life, v), vi) and vii) are satisfied because of the 
presence of the Moon, while conditions v) and viii) are the result of 
the planets all having small eccentricities and mutual inclinations, 
and the Earth being far from a secular orbital resonance and from 
a (secular) spin-orbit resonance. Thus at first glance it appears that 
the Moon plays a key role in supporting life on Earth. However, 
some of the above definitions should be interpreted with care. Mars 
appears to fulfil several of these criteria (moderate rotation speed, 
moderate obliquity, moderate orbital perturbations, moderate 
seasonal temperature variations). In addition, it could have a stable 
spin-axis if the secular frequencies where different. However, the 
slow rotation comes from the fact that it is a planetary embryo 
(Dauphas & Pourmand, 2011), and embryos appear to be less 
adapted for life because they lack the gravitational strength to 
keep a thick atmosphere and oceans against Jeans escape, impact 
erosion or stellar wind pick-up over 4.5 Gyrs. The mass of an 
embryo may be regulated by an 'isolation mass', which is one 
order smaller than an Earth mass at 1 AU. Thus for the remainder 
of this paper we focus on a fully-formed planet with a satellite, 
such as the Earth-Moon system. But is the Earth-Moon system 
unique in its architecture or is the current system a common 
outcome? We aim to answer that question and relate the findings 
to the habitability of the resulting systems. 

It is generally believed that the Moon formed through the 
giant impact of a Mars-sized body with the proto-Earth (e.g. 
Hartmann & Davis, 1975; Cameron & Ward, 1976; Ida et al, 
1997; Kokubo et al., 2000; Canup, 2004). This impact melted 
the impactor and part of the proto-Earth, and resulted in much of 
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the ejecta orbiting the Earth (Benz et al., 1991; Cameron, 1997; 
Canup, 2004). The Moon subsequently accreted from this disc of 
debris in less than a year (Ida et al., 1997). Usually the final disc 
mass was less than 2.5 lunar masses and orbited inside the Earth's 
Roche limit, currently at 2.9 Earth radii (i?e)- The Moon ends up 
orbiting in the Earth's equatorial plane. Tidal evolution over the 
next 4.5 Gyr pushed it towards its current orbit (Goldreich, 1966; 
Touma & Wisdom, 1994). With the aid of numerical simulations 
Elser et al. (2011) conclude that binary planetary systems, such 
as the Earth-Moon system, occur approximately 8% of the time, 
ranging from 2.5% up to 25%. However, their studies did not take 
subsequent tidal evolution into account. 

Goldreich (1966) and Touma & Wisdom (1994) have shown 
that the Earth's current obliquity of 23.5° requires its obliquity to 
have been approximately 10° just after impact. Similarly, Earth's 
rotation period was approximately 5 h after the impact and the 
Moon's inclination with respect to the Earth's equator was about 
10°. The initial spin of the Earth being nearly perpendicular to its 
orbital plane is a fairly rare outcome after a Moon-forming impact: 
it has been shown that the obliquities of the terrestrial planets are 
isotropically distributed (Chambers, 2001; Kokubo & Ida, 2007) 
and their resulting rotation rates are approximately half of their 
maximum values (Kokubo & Genda, 2010). It is likely that the 
planetary rotation period is partially determined by the impact 
that formed the satellite (Lissauer, 2000) and the satellite mass is 
related to the impact parameter of the collision: grazing impacts 
produce heavier satellites and speed up the planet's rotation rate 
much more than head-on collisions (Kokubo et al., 2000). This 
result suggests there is a relation between the planet's rotation 
period and the satellite that has formed. However, as we show 
in Appendix A, the relation is not one-on-one but rather a rough 
estimate. Additionally, an isotropic obliquity distribution favours 
coplanar spins rather than perpendicular ones. What is then the 
evolution of the Earth-Moon system if the Earth had been much 
more oblique just after the Moon-forming impact? 
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Figure 1. This figure shows the regions of the three possible outcomes of 
planetary motion as a function of satellite mass and original planetary obliq- 
uity. Taken from Atobe & Ida (2007). See text for details. 



of Al, where a retrograde-spinning planet will eventually spin 
prograde due to the dominance of solar tides over satellite tides 
(for Al this does not happen). 

From this figure it appears that the outcome for oblique 
planets can be very different from the current Earth-Moon system. 
Atobe & Ida (2007) studied the evolution of these systems until 
their ultimate state. For light satellites, such as the Moon, the time 
scale to reach the final configuration is much longer than the age 
of the Solar System, exceeding the Hubble time for masses below 
half a lunar mass. Thus their results cannot be used to directly 
determine the state of the system at the current epoch i.e. 4.5 Gyr 
after its formation. 



Atobe & Ida (2007) studied the tidal evolution of systems 
consisting of an Earth-like terrestrial planet with a satellite ranging 
from sub-lunar to super-lunar mass on a circular orbit. They varied 
the initial mass of the satellite (m s ) and the initial obliquity of 
the planet (so) but kept the initial semi-major axis of the satellite 
at 3.8 R®, the Earth's rotation period at 5 h and the satellite's 
inclination with respect to the Earth's orbit, i, at 1°. They showed 
that there are essentially three outcomes when these systems 
are evolved towards completion: i) the satellite recedes from 
the planet and approaches (and sometimes reaches) the outer 
co-rotation radius (case A), ii) the satellite first recedes from and 
then approaches the planet resulting in collision (case B), or iii) 
the satellite first recedes from the planet and then approaches it 
but gets caught at the inner co-rotation radius (case C). In general 
terms, case A occurs for low-mass satellites at low prograde or 
high retrograde obliquities. Case B occurs for initial obliquities 
£o € (~ 60° , ~ 120° ) and case C occurs for heavy satellites 
( m s > 0.03 m p , with m p being the planetary mass) at low 
obliquity. We have summarised the results in Fig. [TJ taken from 
Atobe & Ida (2007), where the symbol 70 is eo in our notation. 
The outcome case A is further subdivided according to whether 
the system is still evolving. For case Al prograde planetary spins 
will have the planetary obliquity, e p evolve to e p — > 0° while 
retrograde spins will evolve to e p — > 180°. In the outcome A2 
the system has stopped evolving. Outcome A3 is a particular case 



In this paper we determine the state of fictitious systems 
consisting of an Earth-like terrestrial exoplanet and a satellite 
whose mass ranges from sub-lunar to super-lunar values. We 
investigate the tidal evolution of these systems starting with 
different initial planetary obliquities and satellite masses. We do 
this by performing a series of simulations similar to those done 
by Atobe & Ida (2007) but we cease them after the system has 
reached an age of 4.5 Gyr. We then analyse the states of each 
system and determine which regions of the phase space are the 
most suitable for supporting terrestrial-type life. In addition we 
attempt to determine how likely it is to produce a system similar to 
our own from a range of initial conditions. This paper is structured 
as follows. 

In the next section we lay out the equations of motion con- 
sisting of the star, planet and satellite. These are based on the work 
of Boue & Laskar (2006) and Correia (2009). In section 3 we de- 
scribe our numerical methods and compare the evolution from our 
simulations to earlier results presented in Boue & Laskar (2006) 
for the conservative motion, and Touma & Wisdom (1994) for the 
tidal evolution. In section 4 we present the results of our numeri- 
cal simulations, focusing on the final states of the system and how 
these may change with slightly different initial conditions and pas- 
sages through secular spin-orbit resonances. In the last section we 
present a summary and conclusions. 
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2 THEORY OF PRECESSION AND TIDAL EVOLUTION 

We study the evolution of a hierarchical system consisting of 
a central star with mass M», a terrestrial planet with mass m p 
and a satellite with mass m 3 . The planet and satellite are oblate 
spheroids with zonal harmonics J2 = (C — A)/m p R 2 , where R 
is the equatorial radius of the body under consideration and C and 
A are the principal and secondary moments of inertia. The planet 
and satellite rotate about the axis with maximal inertia, C. From 
now on, we shall use the dimensionless variant C = C/m p R 2 . We 
follow Boue & Laskar (2006) and Correia (2009) in the derivation 
below. 

The conservative motion of the system can be tracked through 
the conservation of the total angular momentum 



L p + L s + H„ + H s = 0, 



(1) 



where a dot denotes a time derivative, H; = CiiriiR^ViSi are 
the spin angular momenta of the planet and satellite and Li = 
rrnma 2 (1 — ei) 1 ^ 2 ki are their orbital angular momenta. Here Vi 
is the spin frequency of either body, cii are the semi-major axes, d 
are their eccentricities and rn are their mean motion. The vectors s; 
are the unit vectors in the direction of the spin angular momentum 
of the planet and the satellite while the vectors ki are the unit vec- 
tors of the orbital angular momenta. In what follows the subscript i 
can refer to either the planet or the satellite. Tidal evolution acts on 
time scales much longer than the orbital period of both the planet 
and the satellite, so equation l[T) is averaged over the mean anoma- 
lies of both bodies and the periapse of the planet-satellite pair. The 
resulting equations for the conservative motion are then given by 
(Boue & Laskar, 2006) 



Hi = — «i COS £j kp X Si — Pi COS 9i k s X Si, 

L s = —7 cos J k p x k s + ^ jij cos 9j k s x s i; 



(2) 



Lr} — 



7 cos I k p x k s + on cos Si k p x Si . 



The equations above assume that the eccentricity of the satellite is 
constant during one full rotation of the line of apses. This does not 
necessarily have to be true, but we shall assume that it is so in this 
simplified model. A study that also accounts for the eccentricity 
evolution of the satellite is left for the future. 

In equations l|2} above we have introduced the angles / = 
arccos(k p • k s ), which is the inclination of the orbit of the satel- 
lite with respect to the planet's orbital plane, Ei = arccos(k p ■ Sj) 
are the obliquities of both bodies with respect to the planet's orbital 
plane and 8i — arccos(k s ■ Si) are the obliquities of both bodies 
with respect to the satellite's orbital plane. In this study we are pri- 
marily interested in the quantities J, e p and partially 6 S . We also 
introduced the precession constants 



a, 



A = 



7 



3GM,miJ 2t Ri 
2a 3 ,(l - e 2)3/2 ' 

3Gm p m s J2iRi 
2a|(l-e!)3/ 2 ' 
3GM,m s a 2 s (2 + 3e 2 ) 



(3) 



8a|(l 



e 2 ) 3 / 2 ■ 

Roughly, the quantities aj are the precession rates of the figures 



of the planet and satellite caused by the torques from the Sun, the 
quantities /3j are the precession rates of the poles of both bodies 
caused by mutual torques between them, and 7 is the precession 
rate of the nodes of the satellite's orbit on the planet's orbit. 

On long time scales the system is subject to tidal forces, most 
notably between the planet and the satellite. In this study we adopt 
the constant time delay model of Mignard (1979, 1980) for its sim- 
plicity, keeping the eccentricity of the planet and satellite constant 
(but not necessarily zero). The tidal forces from the planet and the 
satellite on each other act to reduce their rotation rates. The reduc- 
tion in spin angular momentum is compensated by an increase in 
the orbital angular momentum of the satellite. The averaged equa- 
tions of motion governing the tidal evolution of the planet-satellite 
pair are (Correia, 2009) 



H, = -iA"4X - 6 '°(e s )( Sl + cos6lik s )^ 

2n s X - 8 ' (e s )(l-e 2 ) 1/2 k s ], 
L s = § ^ Ki[XQ 6 ' (e s )(si + cosgjk s )^i 



(4) 



2n s X- 8 '°(e s )(l-^) 



2x1/2, 



where we introduced 



K p = 3k2 P Gm 2 . RpAtpa s 6 , 
K s = 3k2sGm 2 R b a At a a~ e '. 



(5) 



Here k<u is the secular Love number of the planet or satellite, Ati 
is the tidal delay and X c j~ 6 '°(e) and XJ~ 8,0 (e) are Hansen coeffi- 
cients. These are defined from 



exp(im«) = X"' m (e) exp(iZM), 



(6) 



where r is the distance of the planet to the star, v is the true 
anomaly, M is the mean anomaly, I, m and n are integers and 
i 2 = — lis the imaginary unit. The time delay is related to the tidal 
dissipation parameter, Q, via At = (Qv)^ 1 . 

The planet-satellite pair is not isolated and it will be subjected 
to tidal forces from the central star. For simplicity we only incor- 
porated the tidal action from the star on the planet, and thus the 
reduction in the planet's spin angular momentum is compensated 
by an increase in its orbital angular momentum, giving 



H,, 



^N[X 6 '°(e p )(sp + cos£ikp)^ p 
8 -°(e p )(l-e 2 )^ kp ], 



2n p X 



(7) 



with N — 3k2pGM t RpAtpa,p . As the satellite recedes from 
the planet and the rotation of both bodies slows down, their J2 
coefficients decrease because the rotational deformation of their 
figures becomes less severe. The J2 values of the planet and 
satellite are updated according to J2i = ^ksiRi^f /Gnu (Atobe 
& Ida, 2007), where k 3 is the secular Love number, which is 
assumed to be constant. 

Now that we have discussed the basic theory of the tidal and 
conservative motion, we describe our numerical methods below 
and compare them with previous results. 
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Quantity Planet Satellite 





3.005 x 10" 6 M 


0.0025 - 0.05 m p 


R, 


4.2634965 xl0~ 5 AU 


1.161848 xlO" 5 AU 


J', 


5 h (3 h and 7 h) 


18.7 h 




0°-175° 


variable 


e, 


variable 





e 


0.331 


0.394 


k 2 


0.300 


0.02264 




0.95 


1.2 


Qi 


20 


30 




1.0 AU (0.75 AU and 1.25 AU) 


3.8 










i 




0, 1° and 5° 



Table 1. Initial values of various quantities for our numerical simulations. 



3 NUMERICAL METHODS AND TESTS 

In order to determine the evolution of the planet-satellite system 
with various initial conditions, we integrated the equations of 
motion with the aid of computer codes. We wrote three versions. 
The first only integrates the conservative motion without any tidal 
evolution. A second version only includes planet-satellite tides 
while the third version also includes the solar tides. The integration 
was performed using the Bulirsch-Stoer method (Bulirsch & Stoer, 
1966). We checked our code against the results of Touma & Wis- 
dom (1994) and Boue & Laskar (2006), which were all reproduced. 

The input parameters are the masses, obliquities, semi-major 
axes, eccentricities, radii, moments of inertia C, tidal dissipation 
factors Q, Love numbers k 2 , mutual inclination i, and rotation pe- 
riods. The values of J2 were derived from the input values. The unit 
vectors were computed from 

k p = (0,0, 1) T , (8) 
k s = (sin I sin CI, — sin I cos f2, cos/) T , (9) 
s s = (sin e p sin tp, — sin e v cost/), cos e p ) T , (10) 

with s s generated similar to s p . For a close satellite we set fi = ip 
and / ~ e p so that i ~ 0. In Table Q] we have listed all the input 
parameters and their initial values. 

It is well known that integrating the tidal equations backward 
in time using the current tidal dissipation rate causes the Moon to 
fall onto the Earth approximately 1 Gyr ago (Touma & Wisdom, 
1994), so we scaled the final integration time to be 4.5 Gyr to 
obtain the current system state. In most cases we also lowered Q p 
to 10 to hasten the evolution and save CPU time. We made sure 
that this did not affect the final outcome. 

To demonstrate the validity of the code and to show the 
different types of precession without tides, we plot the evolution 
of a fictitious Earth-Moon system for various lunar distances in 
Fig. [2] We have taken the current system and simply decreased 
the lunar semi-major axis, keeping everything else the same. Thus 
the evolution at short lunar semi-major axis is unlikely to be 
representative of the past lunar orbit but serves to illustrate the 
various types of motion of the system. 

The top two panels are for the current system with 
a a = 60 Rq, the middle panels have a 3 = 10 R® and the 



bottom panels are for a 3 = 5 R®. The left panels plot the 
inclination of the lunar orbit with respect to the ecliptic (J) in red, 
the inclination of the lunar orbit with respect to the Earth's equator 
(i) in blue, and the obliquity of the Earth with respect to the ecliptic 
(e p ) in green. In the right panels we plotted the evolution of the 
Moon's nodes on the ecliptic (Q) in red, and the node of the Earth's 
equator on the ecliptic (1/)) in green. We distinguish three types of 
motion, as outlined in Atobe & Ida (2007). On the top panels, the 
Moon's orbit and the Earth's spin precess around k p , which results 
in / and e being nearly constant, and i oscillating between e p — I 
and e p + I with a period equal to the revolution of (~ 18.5 yr). 
On the bottom panels, the Moon's orbit precesses about the 
common angular momentum vector of the spin of the Earth and the 
lunar orbit, s p + k s . Thus, the mutual inclination, i, stays almost 
constant while / and e p oscillate out of phase with a period equal 
to half the precession period of the satellite's orbit on the equator 
of the planet (e.g. Kinoshita & Nakai, 1991). Superimposed on this 
is a long-period precession of k s + s p (bottom-right panel). In the 
middle panels, the precession of both s p and k s is neither around 
their sum nor around k p , but somewhere in between (Atobe & Ida, 
2007). Thus /, i and e p all oscillate with large amplitude. Here we 
find that the precession period of the Earth's axis is about 500 yr; 
Boue & Laskar (2006) found it was 450 yr. 

The tidal evolution of the Earth-Moon system to the present 
is shown in Fig. [3] with initial m s — 0.01223 m^, e p = 7.3° and 
/ = 19° . The values of the other relevant quantities are in TableQ] 
The top-left panel depicts the obliquity of the Earth (green) and 
the inclination of the Moon with respect to the ecliptic (red) versus 
the semi-major axis of the lunar orbit in Earth radii. To reproduce 
the current lunar ecliptic inclination, a high original inclination 
close to the Earth is needed, i ~ 11° (Touma & Wisdom, 1994). 
This appears in contradiction with lunar formation theory from 
an impact (Canup, 2004). A possible solution to this problem is 
if the Moon passed through the evection and eviction resonances 
at 4.6 i?0 and 6 R® (Touma & Wisdom, 1998), or maybe if the 
Earth originally had two moons (Jutzi & Asphaug, 201 1). 

The top-right panel of Fig. [3] depicts the rotation period of 
the Earth vs the semi-major axis of the Moon. The results of both 
the top panels agree with those presented in Touma & Wisdom 
(1994). The bottom-left panel depicts the semi-major axis of the 
lunar orbit with time. It is well known that the current dissipation 
in the Earth is anomalously high (e.g. Williams, 1999), caused 
by a near-resonance between ocean free modes and tidal forcing 
(Webb, 1982). Here we have just scaled the time such that the 
current system is reproduced at 4.5 Gyr. Unlike the tidal model 
with constant phase lag, where Q is independent of the tidal 
forcing frequency, for the Mignard tidal model the semi-major axis 
evolution is not expressible in closed form as a function of time. 
At small semi-major axis a s (t) oc t 15 and for large semi-major 
axis a s (t) oc t , similar to, but not the same as for the constant 
phase lag model where a 3 (t) oc t 2//13 « t 0,153 . In the same figure, 
the bottom-right panel depicts / and e p as functions of time: I 
decreases monotonically to about 5° while e p continues to increase. 

The result of a final test is depicted in Fig. [4] which plots 
the precession frequency of the Earth's spin (black line) and 
the corresponding precession period (grey line) as a function of 
lunar distance. The data points were generated from numerical 
simulations: we tidally evolved the Earth-Moon system up to a 
pre-determined value of the semi-major axis of the Moon, using the 
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Figure 2. Orbital evolution of the Earth-Moon system for several different initial configurations. The top panels depict the current state with the Moon at 
60 K®. The middle panels pertain to the Moon being at 10 H® and the bottom panels have the Moon at 5 H®, On the left column we plot / (red), e p (green) 
and i (blue). The right panels plot Q (red) and ip (green). 



initial conditions of Table Q] with initial e p = 7.3° and I = 19°, 
similar to Touma & Wisdom (1994). Once the Moon had evolved 
to the designated distance, the simulation was stopped and the 
final state was integrated for 100 kyr without tides. The general 
shape and magnitude of the curve agrees with Touma & Wisdom 
(1994) and Boue & Laskar (2006), but the maximum is at a larger 
lunar semi-major axis. We believe this is caused by us changing 
the value of the Earth's J2 while Boue & Laskar (2006) kept it 
constant. 

Now that we have demonstrated the validity of our codes 
through comparison with earlier work, we performed a series of 
simulations similar to Atobe & Ida (2007): we fixed all the in- 
put quantities of the system to their values listed in Table Q] but 
changed the initial obliquity of the planet and the mass of the satel- 
lite. The obliquity was increased from to 175° in steps of 5° 
while the satellite mass ranged from 0.0025 m p to 0.05 m p in 
steps of 0.0025 m p . Each of these planet-satellite systems were in- 
tegrated for 4.5 Gyr and their endstates were recorded. These end- 
states were subsequently integrated for 100 kyr without tides to 
determine the spin precession frequency of the planet and compare 
them with analytical results of Boue & Laskar (2006). These end- 
states should quantitatively reflect the possible configurations of 
the planet-satellite system for different values of the initial satellite 
mass and initial planetary obliquity. While the outcome is similar to 
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Figure 4. Precession frequency, ip (black line) and the period of revolution 
of the Earth's spin axis (grey line) as a function of lunar separation. 



Atobe & Ida (2007), our approach is different because we terminate 
the simulations at a given time rather than determine the final state 
of the system after all tidal evolution. We also checked whether the 
planet's precession frequency will pass through a secular spin-orbit 
resonance with one or more of the inclination eigenfrequencies of 
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Figure 3. Tidal evolution of the Earth-Moon system. This plot should be compared with Touma & Wisdom (1994). The top-left panel plots / (red) and e p 
(green) vs cim ■ The top-right panel shows the Earth's rotation period vs ajvf . The bottom-left panel shows a^j vs time and the bottom-right panel depicts / 
(red) and e p (green) vs time. 



the solar system, s* (Brouwer & van Woerkom, 1950). Such a res- 
onance would induce a large, sudden increase in the planet's obliq- 
uity. We saved the state of the system whenever such a resonance 
was crossed for future analysis. 



4 RESULTS 

In this section we present the results from our numerical simu- 
lations. We have used the third version of our code i.e. the one 
implementing both satellite and solar tides. We show the final 
state of the planet-satellite system at 4.5 Gyr in the form of 
several figures. The results of simulations with different initial 
planetary spin periods and planetary semi-major axis are shown 
further below. Whenever the satellite collided with the planet the 
simulation was stopped and the final obliquity and spin period of 
the planet were recorded. No subsequent solar tidal evolution was 
taken into account. 
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Figure 5. Contour plot of the satellite's semi-major axis (in planetary radii) 
after 4.5 Gyr of tidal evolution. In the central wedge-shaped region the satel- 
lite collides with the planet shortly after its formation (Atobe & Ida, 2007). 



4.1 Endstates 

In Fig. [5] we plot the semi-major axis of the satellite in planetary 
radii after 4.5 Gyr of tidal evolution as a function of the mass of 
the satellite and the initial obliquity of the planet. All subsequent 



figures will be of this type. The colour coding on the right shows 
the scale. The deep purple region in the middle of the figure, which 
flares out with increasing satellite mass, is the region where the 
satellite falls onto the planet (Atobe & Ida, 2007). In all the figures 
the collision region will be this colour. This endstate occurs for eo 
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Figure 6. Colour surface plot of the period ratio of the planet's rotation 
and the satellite orbit. A ratio of 1 means the system is in the double syn- 
chronous state. 

ranging from approximately 60° to 120°. Here extreme seasonal 
variations in insolation and temperature may occur (e.g. Williams 
& Kasting, 1997) and these could be a concern for terrestrial-type 
life. The widest satellite orbits are obtained for low-obliquity 
planets and masses between 0.005r7i p and 0.0125m p , and for 
planets with initial obliquity close to 180°. Smaller satellites do 
not raise a high enough bulge on the planet and thus need more 
time to expand. It should be noted that in all figures, the structure is 
symmetrical around eo = 90° . The size of the satellite orbit affects 
the strength of the tidal bulge on the planet and its pole precession 
frequency. The latter, especially, could be important for habitability 
when this frequency is close to one of the eigenfrequencies of the 
solar system. 

Something different occurs for massive satellites, which is 
clearly seen in Fig. [6] Here we plot the ratio of the rotation period 
of the planet to the orbital period of the satellite. For satellites more 
massive than 0.02 m v the system is in the double synchronous 
state. Atobe & Ida (2007) call this evolution type C. The final pe- 
riod ratio depends steeply on m B and scales as m* 4 (Atobe & Ida, 
2007) and thus the transition from an evolving system to a double 
synchronous one with increasing m s is rather abrupt. Applying 
this to our own system suggests that the habitability of the Earth 
would be completely different if the Moon were just 50% more 
massive. The double synchronous state also explains the decreas- 
ing final satellite semi-major axis with increasing satellite mass: 
the system reaches the double synchronous state and the final orbit 
is closer to the planet because of the lower initial ratio of angu- 
lar momentum in the planet's rotation vs. that in the satellite's orbit. 

The corresponding rotation period of the planet after 4.5 Gyr 
of tidal evolution is shown in Fig. [7] Here we plot log(P r /lhr) 
as a function of satellite mass and initial planet obliquity. Similar 
to Fig. [6] above, there is a very steep increase of the rotation 
period as a function of satellite mass, and the period reaches a 
maximum for a mass of approximately 0.02m p . Here the rotation 
period of the planet is several weeks and the satellite semi-major 
axis sits between 40 and 45 R p . As the satellite mass increases, 
the final rotation period of the planet is reduced because of the 
larger fraction of initial angular momentum residing in the satellite 
orbit compared to the planet's rotation. Even so, in the double 
synchronous state the rotation period of the planet is always longer 
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Figure 7. Colour surface plot of the logarithm of the planet's final rotation 
period in hours. 

than 96 hours, apart from very close to the edge of the central 
wedge. 

As we discussed earlier, these long rotation periods affect 
the climate on the planet and have implications for the origin and 
sustenance of life. In the double synchronous case the rotation 
period of the planet could be a substantial fraction of the relaxation 
time scale and large diurnal temperature variations could occur. 
Furthermore, in the double synchronous state there are no tides 
and thus ocean mixing could stop, thereby reducing possible 
heat transport to the poles. However, the slow rotation rate may 
increase the Hadley cell to reach the poles, which increases the 
heat flow to and the temperature at the polar regions. More studies 
are needed to decide which of these two effects dominates over the 
other. Lastly, the long rotation period reduces the J2 moment and 
causes the precession period of the planet to increase. This issue is 
discussed below. 

Systems with very light satellites (say m s < 0.005 m p ) may 
negatively affect the habitability because the lack of a tidal bulge 
raised by the satellite causes the planet to still rotate rapidly, which 
reduces ocean mixing. This would decrease the heat flow towards 
the poles, likely causing sustained low temperatures there that 
reduce habitability. 

Even though many systems are in the double synchronous 
state after 4.5 Gyr, the real test to measure if the system has fully 
evolved is to determine the final obliquity of the planet. From 
Atobe & Ida (2007) we know that when the double synchronous 
state is reached, the planet's obliquity slowly decreases back to 
if prograde, or increases to 180° if retrograde. Figure [8] depicts 
the obliquity of the planet after 4.5 Gyr of tidal evolution. All 
double synchronous systems are fully evolved and the planet's 
obliquity is at or 180°. The only cases where the obliquity has 
a different value is when the system is not synchronous or when 
the satellite has collided with the planet. For planets whose spin is 
perpendicular to their orbit, the habitable regions are substantially 
reduced from planets with moderate to high obliquities (Spiegel et 
al., 2009). 

One last quantity that determines the stability of the obliquity, 
and potentially the habitability, is the precession frequency of 
the spin pole of the planet. Figure [9] plots the logarithm of the 
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Figure 8. Contour plot of the final obliquity of the planet in degrees. 
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Figure 10. Cumulative distribution of satellite mass generated through a 
giant impact. See text and Appendix A for details. 
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where p So is the probability of the initial obliquity to be £o and 
p m is the probability of the satellite to have a mass m s . The 
numerator sums over all the non-synchronous cases while the 
bottom sum is for all cases. The probability of the initial spin has 
the distribution p eo deo = | sineocfeo (Kokubo & Ida, 2007). We 
find that the total probability that the system is still evolving is 85%. 
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Figure 9. Contour plot of the logarithm precession frequency, in arcsec per 
year, of the planet's spin. Contours for resonances with the eigenfrequencies 
sg and S3 are indicated. 



precession frequency of the spin pole of the planet, in arcsec per 
year, as a function of satellite mass and initial obliquity. In order to 
have a rough idea of whether this system will experience a secular 
spin-orbit resonance, contours corresponding to resonances with 
the nodal eigenfrequencies ss and S3 are indicated. Since S4 ~ S3 
(e.g. Brouwer & van Woerkom, 1950) we did not include it here. 
We realise that in a system different from our own Solar system, 
the eigenfrequencies should not be the same, and thus the contours 
serve an illustrative purpose only. The precession of the planet is 
slowest when the rotation period is the longest, and some of these 
systems could have evolved through the ip = se and tp = S3 and 
tp — S4 secular spin orbit resonances. 

From the above figures it appears there are generally three out- 
comes: i) the system is still evolving, ii) the system is in the double 
synchronous state, or iii) the planet has no satellite (initial plane- 
tary obliquity near 90°). The current Earth-Moon system belongs 
to the first category. We can estimate the probability of being in the 
first state once we know the mass distribution of the satellites from 
giant impacts. We have performed this analysis and detailed it in 
Appendix A and displayed the result in Fig.[l0] The probability for 
the system to be still evolving is the probability of the system to not 
be synchronous, and can be evaluated as 



However, what systems yield a state similar to our own 
i.e. they have 12 < P r < 48 h, e p < 40° (> 140°) and 
0.005 < m s < 0.02 m p ? From examining Fig. [8] we can find 
the initial obliquity that yields a system with the final planetary 
obliquity and spin period in the specified ranges. We can then 
repeat the same procedure as above, which yields a probability 
of 14% for the system to still be evolving, have e p < 40° or 
£ p > 140° and 12 < P r < 48 h. This value is somewhat uncertain 
because there is no unique method to obtain this probability, and 
thus it only serves as a rough estimate. 

There are several effects that we have not taken into account, 
such as the planet's distance from the Sun, the initial rotation period 
of the planet and the effect of the perturbations of the other planets 
on the obliquity of the planet. Each of these will be discussed in 
turn. 



4.2 Secular spin-orbit resonance crossing 

It is unlikely that extrasolar terrestrial planets will exist in isolation, 
and recent data confirms this hypothesis (Lissauer et al., 2012). 
The existence of multiple planets in a given system raises the 
possibility of the system crossing secular-spin orbit resonances. 
The chances of that happening depends on the secular architecture 
of the system. Atobe et al. (2004) studied the effect of the obliquity 
evolution of terrestrial exoplanets that were perturbed by a giant 
planet. They concluded that the terrestrial planet's obliquity vari- 
ations were too large to sustain life if the giant planet was closer 
than about 5 Hill radii from the terrestrial planet. However, their 
study was necessarily limited to a few representative cases of plan- 
etary configurations because a general study is currently unfeasible. 

The only planetary system for which we know the secular 
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Figure 11. Obliquity evolution of the planet when passing through the 
ip = S6 resonance for two different expansion rates of the satellite or- 
bits. The jump in obliquity is consistent with the theory presented in Atobe 
et al. (2004). The large obliquity ranges near the right at the top panel are 
caused by the obliquity crossing the resonances ip = S3 and 1/1 = 54. These 
resonances overlap and the motion is chaotic. 

eigenfrequencies with a decent precision is our own. Ideally we 
would like to test how a planet-satellite system responds to a 
secular spin-orbit resonance crossing in a multitude of systems, 
preferably even in a generalised manner, but that is not possible 
at this stage. Thus in what follows we shall focus on a fictitious 
planet-satellite system crossing a secular spin-orbit resonance in 
our own solar system as a proof of concept, and cautiously use it as 
a possible outcome for other systems. If the mutual inclinations of 
the planets in other systems is small, the outcome presented here 
should be quantitatively similar to what happens in other systems. 

Figure [9] showed that it is possible for the planet-satellite sys- 
tem to cross several secular spin-orbit resonances, mostly ip = sa 
and ip = S3. The effect of such a resonance crossing is to increase 
the planet's obliquity. During a secular spin-orbit resonance the 
obliquity oscillates around its resonant value, e r . While in reso- 
nance, the maximum libration amplitude of the obliquity is given 
by (Atobe et al., 2004) 

I A cos E r I max = V2iVisin2£ r , (12) 

where Ni is the forced inclination on the planet's orbit correspond- 
ing to the frequency Si. If the oscillation range is small, then we 
can use the approximation | Ae r | max = 2^/Ni cot e r . The increase 
in obliquity when crossing a resonance is then ~ 2j Aey |. Typically 
Ni cot Ei — O(10 -3 ) to O(10~ 2 ) depending on the resonant 
obliquity and the forcing, and thus typically Ae r ~ 5° -10°. 

In Fig. QT] we portray the obliquity of the planet vs time as 
it goes through the resonance crossing ip = se. The precession 
and obliquity of the planet were integrated using the method 
described in Brasser & Walsh (2011), and the precession constant 
was artificially enhanced so as to match the precession rate of 
the planet just before it crosses the resonance with sq. During 
the integration the precession constant was linearly decreased in 
order to mimic the regression of the satellite from the planet. In 
the top panel the regression was twice as fast as in the bottom 
panel. One can see that, when crossing the resonance, the obliquity 
jumps by about 10° in a few million years in the bottom panel 
while it completes one resonant oscillation in the top panel. In the 



rightmost part of the top panel the planet hits the resonance ip = S3 
and the obliquity oscillates chaotically with large amplitude. For 
a low-obliquity planet the crossing through the resonance with sg 
could have disastrous consequences if it has polar ice caps such as 
the Earth because the increased insolation at high latitudes would 
partially melt these ice caps, raise sea levels and potentially wipe 
out a substantial portion of coastal life. Thus we argue that for 
terrestrial-type life to develop and be sustained, passage through 
secular spin-orbit resonances should be avoided. 

At present the precession rate of the Earth's spin pole is 
ip = —50.5 " yr . However, a 50% increase in the lunar mass 
or initial obliquity of the Earth would have resulted in the Earth- 
Moon system having passed through the ip = se resonance before 
4.5 Gyr. 



4.3 The effect of initial planetary rotation period 

Recent simulations of terrestrial planet formation with a realis- 
tic accretion scenario demonstrated that these planets have rota- 
tion periods of about twice their minimum value, with P m i n = 
2-k sj I Gm p and a spread of about P m i n (Kokubo & Genda, 
2010). In order to account for the different initial rotation period 
of the planet we performed a series of simulations where we set the 
planet's initial period equal to 3 h and 7 h. Rather than show the 
full results we decided to only show the regions where the system 
is doubly synchronous. The results are plotted in Figsll2landll3l 
In the case the initial rotation period is 3 h the minimum satellite 
mass for which the system is synchronous has now increased be- 
yond 0.03 m p , while it is close to 0.0125 m p for the case where 
the initial rotation period was 7 h. Another interesting feature is 
that the area in the plot where the satellite falls onto the planet also 
depends on the initial rotation period, and it grows as the initial 
period increases. The reason for this behaviour is that the tidal evo- 
lution increases the obliquity of the planet as the satellite recedes 
from the planet. At high obliquity the satellite will eventually turn 
around and fall onto the planet, but this only happens if the incli- 
nation between the spin of the planet and the orbit of the satellite 
i — arccos(s p • k s ) > 90° part of the time. The planet's obliq- 
uity increase is lower for fast-spinning planets, a smaller number 
of cases will experience retrograde motion and the satellite will not 
fall onto the planet. 



4.4 The effect of initial planetary semi-major axis 

The habitable zone around the Sun is thought to reside between 
0.7 AU and 1.3 AU (Williams & Kasting, 1997). The precession 
time of the planet's spin and the satellite's orbit scales with the 
planet's semi-major axis as a p 3 . This strong dependence means 
that at the edges of the habitable zone the strength of the solar 
tides varies by up to a factor of 3 from the value at 1 AU. Thus the 
regression rate of the planet's spin and the satellite's nodes vary 
by the same amount, regressing twice as fast at 0.7 AU compared 
to 1 AU and twice as slow at 1.25 AU compared to 1 AU. Since 
the final state of the system is dominated by the tidal interaction 
between the planet and the satellite rather than by the influence 
of the Sun, we found that the final outcome is very similar to 
what was presented in the figures above. However, the different 
precession rates of the planet's spin and the satellite's nodes are 
noteworthy, in particular in the case where the planet is farther 
from the Sun because the region inside the S6 contour of Fig. IT4lis 
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Figure 12. Colour surface plot of the period ratio of the planet's rotation 
and the satellite orbit. The planet's initial spin period is 3 h. A ratio of 1 
means the system is in the double synchronous state. 
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Figure 13. Colour surface plot of the period ratio of the planet's rotation 
and the satellite orbit. The planet's initial spin period is 7 h. A ratio of 1 
means the system is in the double synchronous state. 

greatly enhanced. For planets very close to the Sun (a < 0.5 AU) 
the tidal forces from the Sun may begin to dominate over those 
of the satellite for lunar mass satellites. However, at these close 
solar distances the Hill sphere of the planet becomes comparable 
to the maximum semi-major axis of the satellite for planets with 
low obliquity and approximately lunar mass satellites, and other 
dynamical effects cannot be ignored. Thus the investigation of 
these cases require a proper N-body treatment. This is beyond the 
scope of the current study. 

We have run a simulation of the tidal evolution of a planet- 
satellite pair where the planet has a semi-major axis 1.25 AU. Most 
of the results are the same as in the previous examples, apart from 
the precession frequency of the planet's spin pole. Since ip oc a~ 3 , 
the precession frequency for an identical system at 1 AU should 
be a factor ~ 2 lower. This is depicted in Fig.[l4]above, where we 
show the precession frequency of the planet's spin after 4.5 Gyr of 
evolution as a function of the satellite mass and initial obliquity. 
This should be compared directly with Fig. [9] As one may see the 
contours for se and S3 occupy a larger area of phase space and 
in some extreme cases resonance passage with S2 is possible. As 
we have demonstrated above, resonance with S3 and S4 causes 



Figure 14. Contour plot of the logarithm precession frequency, in arcsec 
per year, of the planet's spin for systems with semi-major axis 1.25 AU. 
Contours for resonances with the eigenfrequencies Sf>, S3 and S2 are indi- 
cated. 

large-amplitude, long-period oscillations in the obliquity and these 
have drastic consequences for the climate. Thus these resonances 
should be avoided. 



5 CONCLUSIONS 

In the previous section we have presented the results from our 
simplified numerical simulations within a well-defined framework. 
Within this framework the results are robust. We have attempted to 
place these results in the context of habitability of tidally-evolved 
terrestrial planet-satellite systems, although the discussion is some- 
what speculative due to the large uncertainties in the habitable 
conditions. 

We have performed a large sample of numerical simulations of 
the tidal evolution of an Earth-like planet with a satellite. The satel- 
lite's mass and the obliquity of the planet are considered as the two 
free parameters; the remaining ones are modelled after the current 
Earth and Moon. In our simplified model, taken from Goldreich 
(1966), Touma & Wisdom (1994), Atobe & Ida (2007) and Correia 
(2009) the satellite's orbit remains circular. Rather than attempt to 
determine the final end state of the tidal evolution we ended the 
simulations when the system reached an age of 4.5 Gyr. This age 
was determined by requiring that the current Earth-Moon system 
is reproduced from the initial conditions of Atobe & Ida (2007). 
We determined i) which systems are still evolving, ii) which ones 
are in the double synchronous state accompanied by a perpendic- 
ular spin, and iii) which systems have lost their satellites. Systems 
with obliquities 60° < e p < 120° lose their satellites (Atobe & 
Ida, 2007). We find that after 4.5 Gyr only 85% of cases, weighted 
by the satellite mass derived in Appendix A, are still evolving; the 
rest have either lost their satellite or have reached the double syn- 
chronous state. We also discussed habitability as a function of the 
planet's obliquity and rotation period in each end state, in terms 
of diurnal/seasonal temperature and ocean mixing and suggest that 
states ii) and iii) may be less habitable than i). Modelling accretion 
of a satellite from debris formed by a giant impact, we estimated 
the probability for an Earth-mass planet to have the end state simi- 
lar to our Earth-Moon system (12h< P v <48h and £0 < 40° or 
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£o > 140°), which might be favourable for habitability, amounts 
to be only 14%. Elser et al. (201 1) conclude that the probability of 
a terrestrial planet ending up with a heavy satellite ranges from 2% 
to 25% with an average of 13%. Combining these results suggests 
that the probability of ending up with a system such as our own is 
on average only 2%. 
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7 APPENDIX A: SATELLITE MASS DISTRIBUTION 
FROM A GIANT IMPACT 

We predict the distribution of mass of a planetary satellite accreted 
from an impact-generated disc. The disc is the result of a collision 
between a protoplanet and a planetary embryo and the satellite 
will form from this disc. The results below are based on SPH 
simulations of giant impacts (Canup et al. 2001; Canup 2004) 
and N-body simulations of accretion of satellites (Ida et al. 1997; 
Kokubo et al. 2000). The purpose here is to get a rough distribution 
but not to pursue detailed fitting with the simulations. 

Ida et al.(1997) and Kokubo et al.(2000) considered conserva- 
tion of mass and angular momentum during the accretion: 



m d ^ m aC c + m a (13) 
L d ~ L acc + L a ~ m^cc^GmpRp + m s ^jGm p a s , (14) 

where mj and Ld are mass and orbital angular momentum of the 
disc, m acc and L acc are those of the disc materials that are eventu- 
ally accreted to the host planet, m s and L s are those that are finally 
incorporated into the satellite, a s is its semi-major axis, and m p 
and R p are a mass and a physical radius of the planet after the giant 
impact. Here we neglected the disc materials that escape from the 
system. From these equations, 



~ 1 J GrripRp + —^-JGmpCLs 

m d V md) M d 



(15) 



Canup et al. (2001) compiled the data of previous SPH sim- 
ulations of giant impacts with 7 = 7712/(7712 + mi) — 0.3 and 
ANEOS for an equation of state, where mi and 7712 are target and 
projectile masses. They found that md/m p (m p ~ 7711 + 7772) and 
Ld/Lg are given by a function of only Li/L g , where L; is impact 
angular momentum and L g is that of a parabolic grazing impact. 
These are given by 

_ 777 1 7772 , 
Li% — IXiVcscO 



(17) 



being the speed of the impactor at the planet's Hill radius. 

Using the compiled data in Canup et al. (2001), we found that 
the specific angular momentum of the disc is given approximately 
by 

m d /i 

Canup (2004) performed SPH simualtions of giant impacts with 
7 = 0.11-0.15 with the M- ANEOS (Melosh 2000) equation of 
state. The new equation of state results in relatively small L d /m d , 

^--1.2^. (18) 

m d /i 

Note that Canup et al.(2001) and Canup (2004) did not present the 
relations d 1 7fc and j I St . We deduced these relations from their re- 
sults. Substituting eqs. (T7J and J 1 8b into equation l !15t , we obtain 

G/l; ~ ms = i^ bl/3 + (1 " 7)1/3]1/2q - ma ' (19) 

where a ~ 1.5 for the old ANEOS and a ~ 1.2 for the new 
ANEOS. The effects of difference in the equation of state are 
expressed by a slight difference in a. 

From the relations l |17t and dl8b it appears that an impact- 
generated disc would be formed from materials of a projectile that 
do not overlap the target in the line of relative motion. We here 
follow the conventional low-velocity oblique impact scenario (e.g., 
Canup et al. 2001) in order to make the discussion simple, although 
an impact with higher velocity and a steeper angle could also con- 
tribute to formation of a large satellite (Reufer et al. 2012). This 
hypothesis allows us to estimate rrid from simple geometrical argu- 
ments. Here AR = Ri, np — (Ri — R2) expresses the scale of a 
part of the projectile forming the disc and mj can be estimated as 



(AR\ A 1 AR 

p {—) =8 pR *[-KTlh 



AR 



R1+R2 

Rv 



Rt+Rs 



(20) 



where we used that [(Ri + R2)/R P ] is almost constant for 7 
0.1-0.3. We set 

3 



P. 



AR 



md= 8 mp \R 1+ R 2 



where /3 is a constant of O(l). From the definition of AR, 



AR 



Rii 



_ Ri - R 2 _ U_ 
R1 + R2 ~ Ri+ R2 R1+R2 ~ L g S 
where 



(l-7) 1/3 -7 1/3 
(1 - 7) 1 / 3 + 7V 3 ' 



(21) 



(22) 



(23) 



L a = 



777 2 



(Ri + i? 2 )7W = m\/2G(777i + m 2 )(Ri + R2) 



7771 + 7772 

= x/2[ 7 1/3 + (1 - l) l,,i ] 1/2 v^/Gm-p-Rp, 

where Ri is the impact parameter for a two-body encounter, 
R\ and R2 are physical radii of the target and the projectile 
(Ri > R2), jti = 77717772/(7771 + 7772) ~ 7(1 — 7)tt7 p , and we 
assumed that 777 p = mi + 7772 and the internal densities of the 
projectile and the target are the same. In addition from angular 
momentum conservation we have S = yl + / tj| sc with TJoo 



which is ~ 0.14,0.28,0.30 for 7 = 0.3,0.15,0.1, respectively. 
(lB^om eqs. ( 12 1 b and J22b we have 



777£ ^ P_( Li _ £ 
777 p 8 \L g S 



(24) 



If we adopt /3 ~ 1.2, this equation reproduces all the results 
in Canup et al. (2001) (7 = 0.3, the old ANEOS), Canup & 
Asphaug (2001) (7 = 0.1, Tillotson), and Canup (2004) with 
(7 = 0.11-0.15). That is, the scaling relationship in eq. i24\ with 
projectile mass to total mass ratio (7) and a (normalised) impact 
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parameter Li /L g is successful. 

Now we substitute eq. d24t into eq. d!9t to obtain 

m s (3 V2h 1/3 + (1 - 7) 1/3 ] 1/2 " " 1 (U_ _ 



(25) 



™p 8 yJa s /R v - 1 

Since [ 7 1/3 + (1 - 7) 1/3 ] 1/2 ~ 1.2-1.25 for 7 = 0.1-0.3 and Ida 
et al.(1997) shows that a 8 ~ 3.8i? p , eq. l |25b reads as 

3 

(26) 

P \ -^9 L 

with C ~ 0.17-0.25 (a ~ 1.2-1.5). If one were to set 
Li = (2/5)m p RpV, where v is the planet's rotation rate (see 
Section 2), together with the first equation in d 1 6b one has a relation 
between the mass of the satellite and the rotation rate of the planet, 
though it is not a one-on-one relation. For example, we consider 
the case of 7 = 0.13 (£ = 0.31) and ra v — lm®. In this case, 
a collision with Lj = 1.2Lem corresponds to Li/L g = 0.72 



(eq. |16t . Then, eq. \26\ yields m s 
where we used C ~ 0.2. 



0.014m® 



l.lTTlMoon, 



We used equation i25l in a Monte-Carlo method to determine 
the cumulative distribution of satellite masses. We randomly se- 
lected 7 6 (0.05,0.5), Voo was taken from a Maxwellian with 
maximum velocity equal to the planet's escape velocity, and Ri 
was chosen from Ri £ (0.2, 1) in units of R p with R\ being uni- 
form. The resulting distribution is shown in Fig.UOl 
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